
#ifndef _QUICK_YLD_H_
#define _QUICK_YLD_H_

#include "Equations/Momentum_YCoDi.hpp"

namespace Tuna {
  
  template <typename T, int Dim> class Quick_YLD;
  
  template <typename T>
  struct Typeinfo <Quick_YLD<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };
  
  template <typename T>
  struct Typeinfo <Quick_YLD<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };
  
  template <typename T>
  struct Typeinfo <Quick_YLD<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };


template<class Tprec, int Dim>
class Quick_YLD : public Momentum_YCoDi<Quick_YLD<Tprec, Dim> > 
{
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aE;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aW;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aN;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aS;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aF;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aB;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::aP;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::sp;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::dx;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::dy;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::dz;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::bi;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::ei;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::bj;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::ej;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::bk;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::ek;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::dt;    
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::applyBoundaryConditions1D;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::applyBoundaryConditions2D;
    using GeneralEquation<Momentum_YCoDi<Quick_YLD<Tprec, Dim > > >::applyBoundaryConditions3D;
	
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::Gamma;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::Rayleigh;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::phi_0;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::T;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::u;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::v;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::w;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::p;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::dv; 
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::calc_dv_2D;
    using Momentum_YCoDi<Quick_YLD<Tprec, Dim > >::calc_dv_3D;

public:

  typedef Tprec prec_t;
  typedef typename TunaArray<prec_t, Dim >::huge ScalarField;

  Quick_YLD() : Momentum_YCoDi<Quick_YLD<prec_t, Dim> >() { }
  ~Quick_YLD() { }; 

  inline bool calcCoefficients1D() { };
  inline bool calcCoefficients2D();
  inline bool calcCoefficients3D() { };
  inline void printInfo() { std::cout << " Quick_YLD "; }

};


template<class T_number, int Dim>
inline bool Quick_YLD<T_number, Dim>::calcCoefficients2D() 
{
    prec_t dy_dx = Gamma * dy / dx;
    prec_t dx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    prec_t RaGaVol = Rayleigh * Gamma * 0.5 * dx * dy;
    prec_t ce, cem, cep, cw, cwm, cwp;
    prec_t cn, cnm, cnp, cs, csm, csp;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	{
	    ce = ( u(i,j) + u(i,j+1) ) * 0.5 * dy;
	    cw = ( u(i-1,j) + u(i-1,j+1) ) * 0.5 * dy;
	    cn = ( v(i,j) + v(i,j+1) ) * 0.5 * dx;
	    cs = ( v(i,j) + v(i,j-1) ) * 0.5 * dx;

	    if ( ce > 0 ) { cem = 0.0; cep = ce * 0.125; }
	    else {          cem = -ce * 0.125; cep = 0.0; }
	    
	    if ( cw > 0 ) { cwm = 0.0; cwp = cw * 0.125; }
	    else {          cwm = -cw * 0.125; cwp = 0.0; }

	    if ( cn > 0 ) { cnm = 0.0; cnp = cn * 0.125; }
	    else {          cnm = -cn * 0.125; cnp = 0.0; }
		
	    if ( cs > 0 ) { csm = 0.0; csp = cs * 0.125; }
	    else {          csm = -cs * 0.125; csp = 0.0; }

	    aE (i,j) = dy_dx - ce * 0.5 + cep - 2 * cem - cwm;
	    aW (i,j) = dy_dx + cw * 0.5 + 2 * cwp - cwm + cep; 
	    aN (i,j) = dx_dy - cn * 0.5 + cnp - 2 * cnm - csm;
	    aS (i,j) = dx_dy + cs * 0.5 + 2 * csp - csm + cnp; 
	    aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt +
			cem - cwp + cnm - csp;
//			+ (ce - cw) + (cn - cs);
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.

	    sp (i,j) = v(i,j) * dxy_dt - ( p(i,j+1) - p(i,j) ) * dx; 
	    
	    if (i <= ei-1)        { sp (i,j) += cem * v(i+2,j); }
	    else if ( i == ei) { sp (i,j) += cem * v(i+1,j); }

	    if (i >= bi+1)        { sp (i,j) -= cwp * v(i-2,j); }
	    else if ( i == bi) { sp (i,j) -= cwp * v(i-1,j); }

	    if (j <= ej-1)        { sp (i,j) += cnm * v(i,j+2); }
	    else if ( j == ej) { sp (i,j) += cnm * v(i,j+1); }

	    if (j >= bj+1)        { sp (i,j) -= csp * v(i,j-2); }
	    else if ( j == bj) { sp (i,j) -= csp * v(i,j-1); }

	}   
    calc_dv_2D();
    applyBoundaryConditions2D();
    return 0;     
}

} // Tuna namespace


#endif //_QUICK_YLD_H_

















